keim <- c(0.136725664, 
0.111842105,
0.123076923, 
0.083801296,
0.166666667, 
0.134516129, 
0.128, 
0.245783133, 
0.119086022, 
0.142857143, 
0.206741573, 
0.398910082, 
0.523544304, 
0.275870647, 
0.207391304, 
0.12795389, 
0.013268765)

MakeIndex <- function(q1, q0) {
	count.stats <- 1:(q1  + q0)
	index.treat <- combn(count.stats, q1)
	index.cont <- sapply(1:ncol(index.treat), function(j) count.stats[! count.stats %in% index.treat[, j]])	
	index <- cbind(t(index.treat), t(index.cont))
	index
}

MakeRandIndex <- function(q1, q0, n.greps = 9999) {
	q <- q1 + q0
	index <- t(cbind(1:q, replicate(n.greps, sample(1:q, q))))
	index
}

OrdTest <- function(stats, k, index, q1, q0) {
	q <- q1 + q0
	M <- nrow(index)
	r.stats <- matrix(stats[index[1:M, ]], M, q)
	perm.stats <- rowMeans(r.stats[, 1:q1]) - rowMeans(r.stats[, (q1+1):q])
	orig.stat <- perm.stats[1]
	perm.stats <- sort(perm.stats)
	as.numeric(orig.stat > perm.stats[k])
}

IMTest <- function(stats, alpha = .05, q1, q0) {
	q <- q1 + q0
	t.stat <- (mean(stats[1:q1]) - mean(stats[(q1+1):q]))/sqrt(var(stats[1:q1])/q1 + var(stats[(q1+1):q])/q0)
	as.numeric(abs(t.stat) > qt(1-alpha, min(q1, q0) - 1))
}

stats <- c(keim[8:17], keim[1:7])

IMTest(stats, .025, 10, 7)

idx <- MakeIndex(10,7)

# value from Table 1 is .0166 for 5% level test with q1 = 10 and q0 = 7
# two-sided 10% level test, do not reject
OrdTest(stats, ceiling((1-.0166)*dim(idx)[1]), idx, 10, 7)
OrdTest(-stats, ceiling((1-.0166)*dim(idx)[1]), idx, 10, 7)

